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' Abstract. We have undertaken a thorough dynamical investigation of five extrasolar planetary systems using 

extensive numerical experiments on the supercomputer of the Max Planck Institute for Gravitational Physics 
(Albert Einstein Institute). This work was performed as part of the Helmholtz Institute for Supercomputational 
Physics Summer School on "Chaos and Stability in Planetary Systems" (2003). The systems Gl 777 A, HD 72659, 
, Gl 614, 47 Uma and HD 4208 were examined concerning the question of whether they could host terrestrial like 

planets in their habitable zones (=HZ). 

First we investigated the mean motion resonances between fictitious terrestrial planets and the existing gas giants 
in these five extrasolar systems. Then a fine grid of initial conditions for a potential terrestrial planet within the 
Q ' HZ was chosen for each system, from which the stability of orbits was then assessed by direct integrations over a 

^ I time interval of 1 million years. For each of the five systems the 2-dimensional grid of initial conditions contained 

C/3 . 80 eccentricity points for the Jovian planet and up to 160 semimajor axis points for the fictitious planet. The 

computations were carried out using a Lie-series integration method with an adaptive step size control. This 
integration method achieves machine precision accuracy in a highly efficient and robust way, requiring no special 
, - adjustments when the orbits have large eccentricities. 

Itn] ' The stability of orbits was examined with a determination of the Renyi entropy, estimated from recurrence plots, 

and with a more straight forward method based on the maximum eccentricity achieved by the planet over the 1 
million year integration. Additionally, the eccentricity is an indication of the habitability of a terrestrial planet in 
the HZ; any value of e > 0.2 produces a significant temperature difference on a planet's surface between apoapse 
and periapse. The results for possible stable orbits for terrestrial planets in habitable zones for the five systems 
are summarized as follows: for Gl 777 A nearly the entire HZ is stable, for 47 Uma, HD 72659 and HD 4208 
terrestrial planets can survive for a sufficiently long time, while for Gl 614 our results exclude terrestrial planets 
moving in stable orbits within the HZ. 

Studies such as this one are of primary interest to future space missions dedicated to finding habitable terrestrial 
planets in other stellar systems. Assessing the likelihood of other habitable planets, and more generally the 
possibility of other life, is the central question of astrobiology today. Our investigation indicates that, from the 
dynamical point of view, habitable terrestrial planets seem to be quite compatible with many of the currently 
discovered extrasolar systems. 
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1. Introduction 

In this study ^ the extrasolar systems Gl 777 A, HD 72659, 
Gl 614, 47 Uma and HD 4208 were examined on the ques- 
tion of whether they could host additional terrestrial like 
planets in their habitable zones (=HZ). 

Since the discovery of the first extrasolar planetary 
system about 10 years ago (Mayor & Queloz 11995(1 . a 
major point of dynamical investigations has been the de- 
termination of stable regions in extrasolar planetary sys- 
tems, where additional planets on stable orbits could exist. 
Today we have evidence for 105 planetary systems with 
120 planets, where 13 systems have more than one planet. 

Because of the observational methods which are in use, 
our knowledge of extrasolar planets is highly biased: al- 
most all planets were detected via ground based measure- 
ments of the central star's radial velocity. This method 
favors the detection of very large (Jupiter-sized) planets 
that move very close to the main star. Several space mis- 
sions are planned (i.e. COROT, KEPLER, TPF, DARWIN) 
that will be able to find planets with a much smaller mass 
(comparable to that of the Earth) by using other detection 
methods. The first results of these missions will only be 
available after more than a decade, however, theoretical 
studies can currently make predictions on the existence of 
stable terrestial planets within the known systems. This 
serves to establish testable theories of planet formation 
and evolution, as well as aiding the future searches for 
such habitable terrestial planets. 

When one looks at the distribution of semimajor-axes 
of all known exosolar planets, one can see that they are 
located from ^ 0.02 AU to ^ 6.5 AU. These distances 
cover the so-called "habitable zone" (=HZ)^ of a typi- 
cal main sequence star, ranging from ~ 0.7 AU to ~ 1.5 
AU (Kasting et al. ll993|l . Some of these large planets move 
in the region around 1 AU that overlaps with the HZ, and 
leave no room for possible terrestrial planets orbiting the 
host star within the HZ - only a large satellite (like Titan) 
or a Trojan like planet could have all the properties which 
are necessary for orbital stability and also the development 
of a stable atmosphere in such cases. When the large plan- 
ets move well outside (or inside) the HZ they will perturb 
any additional undetected planet that may move within 
this zone. 

For extrasolar systems with only one planet the pos- 
sible existence of additional planets has been investigated 
(e.g. Dvorak et al. l2003bl Pal & Sandor lM?^ . In an exten- 
sive paper (Menou & Tabachnik l2003|l the stability of pos- 

Send offprint requests to: R. Dvorak, e-mail: 
dvorak@astro.univie.ac.at 

^ The results presented here were derived during the 3'^'* 
Helmholtz Summerschool at the Helmholtz Institute for 
Supercomputational Physics in Potsdam: " Chaos and Stability 
in Planetary Systems" from September 1*^' to September 26*'' 
2003 

^ The region were a planet fulfills all the requirements of 
being able to develop life and where liquid water can be present 
on the surface of the planet. 



sible terrestrial planets was studied for all extrasolar sys- 
tems and a classification according to their stability was 
established, but there the role of resonances was underes- 
timated. Resonances - mean-motion resonances and secu- 
lar resonances - can either stabilize the motion of a planet 
and protect it from close encounters and collisions or, in 
other cases, intensify the gravitational perturbations and 
therefore destroy the stability of the orbits. Resonances 
between the orbits of the observed large planets have been 
studied (Callegari et al. lMI^ HadjidemetriouEMI Erdi 
& Pal 20(S|), but a detailed study of the perturbations on 
possible additional planets due to resonances is still lack- 
ing. Thus, any study of the dynamics inside the HZ of an 
extrasolar system has also to include an investigation of 
the resonances. 

Determinations of the dynamical stability of mul- 
tiple planetary systems (Kiseleva et al. l2UU2al 12002 bl 
Beauge & Michtchenko I2003|l and of planets in double 
stars (Pilat-Lohinger & Dvorak r2()()2al Pilat-Lohinger et 
al. l2002b[IM7!^ Dvorak et al. l2()()3ai Lammer et al. lM)^ 
have been done recently. 

In this article we investigate the dynamics inside 
the habitable zone of five dilferent extrasolar systems: 
Gl 777 A, 47 Uma, HD 72659, Gl 614 and HD 4208. Except 
47 Uma, all these systems consist of the central star and 
a Jovian planet. According to the spectral type of the 
main star, the habitable zone was determined; a detailed 
analysis of the mean motion resonances inside the HZ was 
carried out and the full width of the HZ was investigated 
according to the stability of possible additional terrestrial 
planets. This article is organized as follows: Section [3 de- 
scribes the five extrasolar systems for which we performed 
a dynamical study for additional planets, in section |3| we 
introduce the dynamical models and the methods used to 
analyze the time series data resulting from our numerical 
experiments. Scction^lshows how we dealt with the terres- 
trial planet orbits in the major mean-motion resonances 
of the Jovian planet (s). Subsequent sections focus on each 
of the 5 systems individually, detailing our investigations 
and results. We conclude with a summary of these results 
and their implications, including comparison to earlier, 
less detailed studies by Menou & Tabachnik f2003"). The 
appendicies include concise descriptions of the analytical 
and numerical methods used, and may serve as a useful 
reference to the reader. 

2. Description of the Extrasolar Systems 

Tabled shows the main characteristics of the five extraso- 
lar planetary systems which we investigated in this study. 
Except the system 47 Uma (which hosts two planets), all 
consist of a giant planet moving outside the HZ of the 
star; this is a situation similar to our Solar System (=SS) 
when we take into account only Jupiter. Thus, an appro- 
priate dynamical model is the restricted three body prob- 
lem (star -I- planet -I- massless terrestrial planet). A dif- 
ferent dynamical model was taken for the system 47 Uma, 
which looks quite similar to our SS when only Jupiter 
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Table 1. Parameters of the investigated systems - 
obtained from http://www.obspm.fr/planets (1. Sep. 
2003) 



Name 


Mass 


Semimajor 


Eccentricity 






axis [AU] 




Gl 777 A (G6IV) 


0.90 Mq 






Ol 777 A b 

vjri Iff ij 


1 33 Mt 


4.8 


48 + 2 


HD 72659 (GOV) 


0.95 Mq 






HD 72659 b 


2.55 Mj 


3.24 


0.18 


Gl 614 (KOV) 


1.00 Mq 






Gl 614 b 


4.89 Mj 


2.85 


0.38 


HD 4208 (G5V) 


0.93 M© 






HD 4208 b 


0.80 Mj 


1.67 


0.05 


47 Uma (GOV) 


1.03 Mq 






47 Uma b 


2.54 Mj 


2.09 


0.061 ±0.014 


47 Uma c 


0.76 Mj 


3.73 


0.1 ±0.1 



and Saturn are considered. Both 47 Uma and the SS have 
the more massive planet orbiting closer to the star with 
the second planet at approximately twice the distance. 
Furthermore, the 47 Uma planets have roughly a 3:1 mass 
ratio and exhibit low eccentricity orbits very much like 
Jupiter and Saturn in our own SS. The restricted four 
body problem, known to be a good model for an asteroid 
in the SS (taken as a massless body), additionally con- 
siders the secular resonances acting on the motion of a 
massless test body and will also serve as a good model for 
terrestrial planets within 47 Uma. 

In a wider sense we can say that all the central stars are 
solar like stars with masses slightly smaller or larger than 
the Sun. The spectral types for these main sequence stars 
(GO to KO) allow HZs in the range between 0.5 AU and 
1.5 AU. With the exception of HD 4208 and 47 Uma all 
the Jovian planets' orbits have significantly larger eccen- 
tricities than that of Jupiter with sometimes very large 
uncertainties. On the other hand, we expect the deter- 
mined semimajor axes to be rather precise. Our parameter 
study will take into account this discrepancy in the pre- 
cision of the measured orbital elements of the cxoplanets 
under study and will be detailed in the following section. 
The minimum masses of the Jovian planets lie exactly in 
the range we expect for gas giants. With regard to the un- 
known inclinations of these extrasolar systems, we always 
took the masses of table as correct, assuming that the 
orbital plane of the system is seen edge-on. It should be 
kept in mind that the true masses of these planets could 
be significantly larger. Statistically we can say that in 5 
of 6 cases the inclination will change the mass by a fac- 
tor smaller than four compared to the one given in the 
tables. Recent computations have shown that with these 
larger masses the strenghts of the perturbations will not 
significantly change (Sandor et al. i2004(l . 



3. Simulation Method and Stability Analysis 

The availability of a supercomputer with 128 processors'^ 
for this investigation enabled the direct computation of 
orbits to assess stability. Futhermore, the extent of the 
computational resources favoured the use of a very precise 
numerical integration scheme, the Lie-integration method, 
which is free from numerical difhculties experienced by 
other (lower order) techniques, particularly in the case of 
highly eccentric orbits. The Lie-integration method uses 
an adaptive stepsize and is quite precise and fast, as has 
been shown in many comparative test computations with 
other integrators such as Runge-Kutta, Bulirsch-Stoer or 
symplectic integrators. Although, symplectic integrators 
are very effective when eccentricities remain small, the 
Lie integrator is the better choice in studies such as this 
one, where very large eccentricity orbits are explored. We 
detail this integration method in Appendix^ and further 
details can be found in Hanslmeier & Dvorak 19811 as well 
as Lichtenegger 1198^ 

The uncertainties in the orbital elements derived from 
observations make it sufficient to consider only Newtonian 
forces in a dynamical model of two massive bodies (the 
central star and the discovered planet) and a massless test 
body, representing a potential terrestrial planet. As men- 
tioned previously, this system corresponds to the three 
dimensional restricted 3-body problem. For the 47 Uma 
system where two Jovian planets are known, we followed 
the dynamics of a restricted 4-body problem with a cen- 
tral star, both discovered Jovian planets, and a massless 
test body representing a potential terrestrial planet. This 
increased the computing time considerably for 47 Uma 
compared to the other restricted 3-body systems. 

For the integrations we use a fine 2-dimensional grid 
of initial conditions; for the terrestrial planets varying the 
distance to the central star, exploring the HZ, and for the 
(inner - 47 Uma) Jovian planet varying the eccentricity 
of its orbit. We vary the eccentricity of the Jovian planet 
because the values derived from observations have large 
uncertainties and yet this is an important parameter 
for stability of other planets within the system. Table |21 
outlines the ranges for our initial conditions. We always 
started the initial orbits of the terrestrial planets as be- 
ing circular; they turned out to deform quite quickly, ex- 
ploring a range of more eccentric orbits so that varying 
their initial eccentricities would seem to be extraneous. 
Besides the eccentricity of the fictitious planet's orbit, 
guided by the formation scenario for rocky terrestrial plan- 
ets (Richardson et al. I2000|l the inclination was set to al- 
most zero. As a final detail we mention that the integration 
of the terrestrial planet's orbit was always started on the 

^ The PEYOTE cluster at the Max Planck Institute for 
Gravitational Physics (Albert Einstein Institute): www.aei- 
potsdam.mpg.de/facilities/public/computers.html 

* In most cases they can only be determined from the shape 
of the periodic stellar velocity curve, while M sin i and the 
semimajor axis come from the zero'th order characteristics, 
the amplitude and period of this curve respectively. 
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Table 2. Initial semimajor axes of the fictitious terrestrial 
planets ax and the eccentricities of the jovian planets ec 



System 




[AU] 






Gl 777 A 


0.5 


- 1.3 


0.4 - 


-0.5 


47 Uma 


0.5 


- 1.5 


0.0 - 


- 0.12 


HD 72659 


0.4 


- 1.6 


0.08 


- 0.3 


Gl 614 


0.5 


- 1.5 


0.3 - 


-0.5 


HD 4208 


0.55 


- 1.4 


0.0 


-0.2 



connecting line between the star and the Jovian planet. 
An integration time of 1 million years was chosen, which 
is long enough to unveil the stability character of the plan- 
ets for the 5 systems of this study. The default parameter 
grid of 80 by 80 which we used for the 5 systems leads to 
a grand total integration time of several 10^" years, when 
we also include test computations and the fact that some 
of the orbits were integrated for an extended 10^ years. 

For the analysis of stability we used - on one hand 
- a straightforward check based on the eccentricities. For 
this we examined the behavior of the eccentricity of the 
terrestrial planets along their orbit and used the largest 
value as a stability criterion; in the following we call it 
the maximum eccentricity method (=MEM). This simple 
check has already been used in other studies of this kind 
and was found to be quite a good indicator of the stability 
character of an orbit (Dvorak et al. I2003a|l . An orbit was 
deemed unstable when the eccentricity exceeded a value 
of e = 0.5, after which we stopped further computation. 
In all our former studies this stability limit turned out to 
be an appropriate tool because all the terrestrial planet 
orbits with e = 0.5 turned out to suffer, sooner or later, 
from a close encounter with the large planet, causing the 
terrestrial planet to escape (Dvorak et al. ■2003al . For the 
habitability of a planet we also used an additional crite- 
rion based directly on the eccentricity of the orbit within 
the HZ. This was done in order to take into account the 
variations in the "solar" insolation on the surface of the 
terrestrial planet. To good approximation (Lammer 2004; 
private communication), requiring e < 0.2 is sufficient to 
keep this variation in insolation small enough during an 
orbit. 

On the other hand we computed the Renyi entropy, 
a measure which is often used in nonlinear dynamics to 
determine how predictable an orbit is. We estimated the 
Renyi entropy by means of Recurrence Plots (RPs), a tool 
of data analysis that has found numerous applications in 
many different fields in the last years (Webber et al. 119941 
Marwan et al. 12002 Thiel et aL l^nHSjl . RPs were initially 
introduced to simply visualize the behaviour of trajecto- 
ries in phase space. The distances of every pair of points of 
the system's trajectory are represented in a 2-dimensional 
binary matrix. In this way RPs yield different patterns 
depending on the system's character. One can introduce 
different measures that quantify the obtained structures. 
It was shown in Thiel et al. (^2003,1 . that it is also possible 



to estimate the Renyi entropy based on the distribution 
of diagonal lines obtained in the RP. The details of this 
second approach are covered in Appendix IXI 

4. The Stability within Resonances 

For the investigation of the resonances, we choose ini- 
tial conditions placed in the most relevant mean-motion 
resonances (=MMRs) of the fictitious planet with the 
Jovian planet inside but also outside the HZ. These res- 
onances were checked for stability in 8 different posi- 
tions of the terrestrial planet (corresponding to M = 
0°, 45°, 90°, 135°, 180°, 225°, 270°, 315°). Additionally the 
computations were carried out with the Jovian planet ini- 
tially placed at the apoastron and periastron. For a de- 
tailed list of the resonant positions that were investigated 
for every system see table 13 




2''0° 



Fig. 1. Schematic view of the stability of orbits in the res- 
onances of I''' and 2"^^ order in HD 4208. Full (empty) cir- 
cles stand for stable (unstable) orbits in apoastron and pe- 
riastron position. When the stability is different we mark 
the apoastron by a triangle, the periastron with a square. 



As an example we discuss the results of the investiga- 
tion of the MMRs for the system HD 4208. We studied 
the following mean motion resonances up to the fourth 
order: 2:1, 3:2 and 4:3 (first order); 3:1, 5:3, 7:5 (second 
order); 4:1, 5:2, 7:4, 8:5, 10:7 and 11:8 (third order); 5:1, 
7:3, 11:7, 13:9, 15:11 (fourth order). As shown in the pic- 
ture for the first set of resonances (1*** and 2""^ order) the 
orbits close to the central star, which move well inside 
the HZ, are aU stable (figure P). For the MMRs close to 
the Jovian planet we can see a preference for stable orbits 
for the initial conditions M = 0°, and 180°. For the 3''^ 
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Table 3. Stability of orbits in mean motion resonances. The numbers give the stable orbits according to the 8 different 
initial conditions 





Gl 777 A 


47 Uma" 


HD 72659 


Gl 614'' 


HD 4208'= 




P 


A 


Mode I 


Mode II 


P 


A 


P 


A 


P 


A 


5:1 


1 


n 
u 


Q 
O 





Q 
O 


Q 
O 


n 
u 


n 
u 


Q 
O 


Q 
O 


4:1 








3 





8 


8 








8 


8 


3:1 














4 


3 








8 


7 


5:2 








6 


7 


2 


2 








3 


2 


7:3 








1 


1 


1 











8 


8 


2:1 








1 


2 


3 


3 








8 


8 


5:3 












2 


1 






7 


5 


3:2 


2 


1 









2 








7 


8 


4:3 








4 













3 


1 


Sum [%] 4.7 1.6 
Total Sum [%] 3.1 


28.6 69.6 
34.9 


42.2 43.8 
43.0 


0.0 0.0 
0.0 


83.3 76.4 
79.9 



" Note that in the case of the 47 Uma system, where two Jovian planets are known, we did not use peri- and apoastron 
position as initial conditions, but 2 different modes corresponding to an aligned or anti-aligned configuration of the two major 
bodies. 

' Besides the 7 given resonances, we calculated the motion inside the 7:2, 9:2 and 8:3 resonance - again, we only found unstable 
motion. 

'= For this system, we calculated all resonances up to the 4*'^ order (see section 4); with the exceptions of the 15:11 and the 
13:9 MMRs the other resonant positions showed the same amount of predominantly stable motion. 




- 8:5 

- 7;i 



Fig. 2. Schematic view of the stabiUty of orbits in the 
resonances of 3"^ order in HD 4208. Description Uke in 
figure n 



90° 




-15:11 
-13:9 
-11:7 



225° 



270° 



Fig. 3. Schematic view of the stabiUty of orbits in the 
resonances of 4*'' order in HD 4208. Description Uke in 
figure n 



order resonance (figure ^ the picture is very inhomoge- 
neous; what we can see is, that for the Jovian planet in 
the apoastron position the orbits are stable even for the 
resonances close to the giant planet. The 4**^ order reso- 
nances are not destabilizing an orbit as we can see from 
figureEl most of them are stable! The percentage of stable 
orbits in resonances is very large for HD 4208. 



Details of the results from the investigation of reso- 
nances for all systems can be found in table O For the 
resonances acting in three of the systems wc can say that 
far from the perturbing planet almost all of them are sta- 
ble; closer to the perturbing planet, they are more and 
more unstable (47 Uma, HD 72659 and HD 4208). Two 
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systems are very much dominated by unstable motion in 
resonances: Gl 777 A and Gl 614. 

5. Gl 777 A 

The first discovery of a planet in Gl 777 A (=HD 190360) 
was reported by Naef et al. H2U(J3|I from the Geneva group 
of observers. This extrasolar planetary system is a wide 
binary with a very large separation (3000 AU); for our 
dynamical investigations of motions close to one star there 
was no need to take into account the perturbations of 
the very far companion. The central star is of spectral 
type G6 IV with 0.9Mq and has a planet of minimum 
mass 1.33 Mj with a semimajor-axis of 4.8 AU. Because 
of the large eccentricity (e=0.48) the possible region of 
motion for additional planets is confined to a < 2.4 AU (= 
periastron). Nevertheless, to have a global stability picture 
of possible additional planets, we investigated the stability 
in the region of the MMR from the 4:3 to the 5:1 resonance 
located at a=1.64 AU. From table|2|one can see that only 
a few percent of the orbits started in the MMRs are stable. 



F 4 ; F 
11 1 z 


1 ! 

i 1 


e 


1- 


GI777AMI 

1.33 


































1 


1 








1 J 

J i l~ 





Fig. 4. Main characteristics of the extrasolar system 
Gl 777 A. The light grey region shows approximately the 
position of the HZ; the dark grey bar indicates how closely 
the planet approaches the central star in its orbit. 



The interesting region of habitability (see figure ^ , 
where planets could have temperature conditions to al- 
low liquid water on the surface, corresponds roughly to 
0.7 < a < 1.3 AU, where we ignore the eccentricity of 
the terrestrial planet. We have started our computations 
in a larger region (0.5 < a < 1.3) with a grid spacing of 
Aa — 0.01 AU and changed also the eccentricity of the 
known planet between 0.4 < e < 0.5 with a gridsize of 
Ae = 0.01. The results of the two methods of analysis of 
the orbital behavior are shown in figures El and El In the 
first plot we show the results of the MEM, where two fea- 
tures are immediately visible: 1) strong vertical lines due 
to high order resonances, and 2) unstable orbits due to 
high eccentricity and high semimajor axes values (red or 
yellow colors). The latter feature is easy to understand 
because closer to the existing planet the perturbations 
are larger. The two methods complement each other in 
the information they convey; the MEM tells us about the 
variable distance to the central star and consequently it 
is a direct measure of the differential energy fiux (inso- 
lation) on the planet. We can therefore determine where 
the variation of this distance does not exceed 50 percent. 



corresponding to an eccentricity of e = 0.2. The Renyi en- 
tropy is a more sensitive probe of the dynamical character 
of the orbit, giving us a measure of the degree of chaotic- 
ity. In particular high order resonance features are made 
very clear using this second method, and we can even see 
the resonances acting when the eccentricity of the planet 
is as low as e = 0.4 (the bottom of figure 0). 

As a result for habitability of a terrestrial planet inside 
the orbit of the Jovian planet, we find that for the system 
Gl 777 A there is quite a good chance that planets will last 
long enough in the HZ to acquire the necessary conditions 
for life in the region with a < 1 AU. 




0.46- 

I 

0.44- 
0.42- 

0.4 -f'^ , , , , , , , 

0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 

(AU) 




eccentricity 

Fig. 5. Initial condition diagram for fictitious planets in 
the system Gl 777 A: initial semimajor axes of the planet 
versus the eccentricity of the Jovian planet. The maximum 
eccentricity of an orbit during its dynamical evolution is 
marked with different colors. 

6. 47 Uma 

The first planet in 47 Uma was discovered by Butler & 
Marcy H1996|l . and the second, the outer one, was reported 
in a paper by Fischer et al. H2002|l . Former studies inves- 
tigated the zones inside the orbit of the more massive in- 
ner planet (Jones & Sleep (21151 Cuntz et al. lMTHjl which 
may host terrestrial planets. The Jovian planet has an 
almost circular orbit with very small errors in the eccen- 
tricity. Because of the errors for the eccentricity of the 
inner Jovian planet's orbit, we varied its eccentricity from 
< e < 0.12 for our computations. It turned out that 
the system became unstable when Cinncr > 0.12 and con- 
sequently also no inner terrestrial planet would survive. 
Almost all the MMRs of the fictitious planet with the 
inner planet are in the HZ (figured). The respective com- 
putations, summarized in table show that almost 1/3 
of planets located in these resonances would survive with 
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0.5- 



0.48 



0.46-J 



0.44- 



0.42-1 



0.4- 
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entropy plot does not show any features other than the 
ones we can see from the MEM. 
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Fig. 7. Main characteristics of the extrasolar system 
47 Uma. Specifications hke in figure 01 
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Fig. 6. Initial condition diagram for fictitious planets in 
the system Gl 777 A: initial semimajor axes of the planet 
versus the eccentricity of the Jovian planet. The value 
of the entropy (=entropy plot) of an orbit is marked in 
different colors. 



moderate eccentricities, which can be explained by the 
small eccentricity of the Jovian planet. The situation is 
similar to Jupiter and Saturn but with a scaling factor of 
approximately 0.4 in a. Also in our Solar System, aster- 
oids are unstable in the outer main belt at and beyond 
the 2:1 resonance, although some islands of stability ex- 
ist, like the Hilda asteroids in the 3:2 MMR. The vertical 
line in figure |H1 located near 1 AU, indicates unstable or- 
bits that mark the 3:1 resonance. The thin vertical line of 
large eccentricities close to a=0.82 AU corresponds to the 
4:1 mean motion resonance of the terrestrial planet with 
the inner large planet. We can observe a broad curved 
"unstable" line between 0.8 AU and 0.92 AU, which is 
caused by secular resonances (=SR) of the perihelia; this 
was confirmed by test computations, where we omitted 
the outer planet. It is remarkable that these SR act very 
strong when the initial eccentricity of the inner planet is 
close to zero, then it becomes weaker and disappears for 
0.08 < e < 0.1, but it is again visible for e > 0.1. It 
is interesting to note how the period of the periastron of 
the giant planet changes (4000yrs < lo < 20000yrs) with 
increasing eccentricity. Between the SR and the 3:1 reso- 
nances a small region of orbits with small enough eccen- 
tricities allowing habitability (e < 0.2) may survive. On 
the outer part of the HZ the very strong 2:1 resonance 
(broad line close to 1.3 AU) limits the region of stable or- 
bits, but still leaves a large region of stable orbits between 
the 3:1 and the 2:1 resonance. This is in good agreement 
with the results of the previously mentioned studies (e.g. 
Cuntz et al. I2003|l . which did not cover such a large vol- 
ume of phase space of possible motions in the HZ. The 
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Fig. 8. Results of the MEM for 47 Uma. 



7. HD 72659 

The G05 star HD 72659 was found to have a compan- 
ion from the Keck Precision Doppler survey (Butler et 
al. .2002^) . The Jovian planet (2.55 Mj) has an orbit with 
a semimajor axis a — 3.24 AU and an eccentricity of 
e = 0.18. The MMRs are located from 2.47 AU (3:2) to 
1.1081 AU (5:1); the 5:1, 4:1 and 3:1 are weU inside the 
periastron position of 2.657 AU and lie in the HZ (around 
1 AU, see figure EJ- 

The resonances turned out to be stable in more than 
40% of the orbits; especially the high order resonances 
close to the HZ are stable in both initial conditions (peri- 
astron and apoastron position). As a consequence we ex- 
pected that these planetary systems may host additional 
terrestrial planets in stable orbits. 

Because of the uncertainties of the observed Jovian 
planet's eccentricity we varied it from 0.08 to 0.30 with a 
stepsize of 0.22/80 = 0.00275 and chose the initial semima- 
jor axis of the fictitious terrestrial planet to satisfy 0.4 AU 
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Fig. 9. Main characteristics of the extrasolar planetary 
system HD 72659. Specifications like in figure 0] 



< a < 1.2 AU. The results are shown in figure [TUl (MEMl 
and figure [TTl fentropv plot). We can identify quite well in 
these plots the resonances up to the 7:1 resonance (only 
in the entropy plot). Again one can see that the dynamics 
of a single orbit can be determined quite well with this 
method; it does not only confirm what is depicted in fig- 
ure^! it also shows many more details especially for the 
motions in resonances. On the contrary the MEM is the 
appropriate tool for determining the eccentricity, which is 
- together with the semimajor axes ~ the crucial param- 
eter for our research of determining planets in habitable 
zones. 
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Fig. 10. Results of the MEM for HD 72659. 



Globally we can see a quite stable HZ in this extrasolar 
system which allows planets on orbits with small eccentric- 
ities. The strong unstable line close to 1.6 AU corresponds 
to the 3:1 resonance, while the other resonances, although 
giving rise to large perturbations in the eccentricities, are 
confined to the center of the resonance. In figure ^| we 
can also see that for the most probable eccentricity value 
of the Jovian planet (e — 0.18) all orbits up to the 3:1 
resonance are stable with low eccentricities (e < 0.2) and 
as a consequence the HZ could be populated with a ter- 
restrial planet (or even more planets depending on their 
masses). 
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Fig. 11. Entropy plot for HD 72659. 



8. Gl 614 

A planet around the KOV star Gl 614 (=14 Herculis) was 
discovered by Naef et al (2004) which is quite massive 
(4.89 Mj) and orbits the central star with a semimajor 
axis a = 2.85 AU and an eccentricity of e = 0.38. The 
computation of orbits in MMRs showed that none of them 
are stable up to the 6:1 resonance in apoastron and peri- 
astron of the existing planet (for details see table O; this 
is not surprising when one looks at the large mass and the 
semimajor axis combined with the large eccentricity: the 
apoastron of the planet is 1.76 AU. The characteristics of 
this systems are shown in figure [T^ 
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Fig. 12. Main characteristics of the extrasolar planetary 
system Gl 614. Specifications like in figure^ 



Looking at the results from the MEM (figure [TS|) one 
can see that for the most probable value of the eccentricity 
of the Jovian planet only the orbits with a < 0.7 stay in 
the stability and habitability (e < 0.2) limits within the 
HZ. 



9. HD 4208 

Using the data of the Keck Precision Velocity Survey, Vogt 
et al. H2002|l discovered that the system HD 4208 contains 
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Fig. 13. Results of the MEM for Gl 614. This plot shows 
that for this system it is very unlikely for a terrestrial 
planet to survive in a stable orbit which also permits liquid 
water to exist on its surface. 



a planet with a slightly smaller mass than Jupiter on an 
almost perfectly circular orbit at 1.67 AU. The G5V star 
has a mass almost like our Sun. Our simulations of the 
dynamical evolution of terrestrial planets were carried out 
in the same way as for the other four systems but with 
special emphasis on the resonances in this system (see 
section 3). Unfortunately, we made an error in the mass 
of the central star (we took it to be 14 % smaller than the 
one determined by the observers) for the computations in 
the HZ (not for the stability in resonances!). The habitable 
zone will be shifted outwards only a few percent because 
of this mass error; motions in the resonances (in the HZ) 
will be, for all practical purposes, unaffected. 

The MMRs of first, second, third and forth order were 
investigated for this system. For a total of 18 resonances 
in two different positions of the existing planet (apoastron 
and periastron) for 8 different values of the mean anomaly 
of the fictitious planets we computed the stability (for a 
discussion see section O. 
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Fig. 14. Main characteristics of the extrasolar planetary 
system HD 4208. Specifications like in figure 0] 



We expected unstable motion for this system from a = 
1.40 AU on, a limit which we computed using estimations 



derived by Wisdom H1980|l and Duncan et al. I|1989(l for 
the onset of global chaos. We took the mean value between 
Wisdom's and Duncan's estimate (namely Aa = 1.27/i^/^ 
- where ^ is the mass ratio of the primaries - which leads 
to a Aa = 0.28 for this system). In the two plots CfigureslT^l 
and I16|l one can see that in fact the orbits are unstable 
from a = 1.3 AU on. 
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Fig. 15. Resuhs from the MEM for HD 4208. 
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Fig. 16. Entropy plot for HD 4208. 

From the 2:1 resonance inward, for all initial conditions 
of the Jovian planet, we can expect stable orbits with low 
enough eccentricities (e < 0.2) to allow habitable terres- 
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trial planets with liquid water, the only exception being 
the resonant orbits (the vertical lines in both figures). 

10. Conclusions 

We have undertaken a dynamical study of five extrasolar 
planetary systems using extensive numerical experiments 
to answer the question of whether they could host terres- 
trial like planets in habitable zones. For the single-planet 
systems we used the elliptic restricted three body prob- 
lem, and for the system 47 Uma the mutual perturbations 
of the two Jovian planets were taken into account and 
therefore the restricted four-body problem served as a dy- 
namical model. Because of the dependence of the stability 
of an orbit in a mean motion resonance on angular posi- 
tion we have selected 8 different positions for the mean 
anomaly of the planet with a step of 45° in the apoastron 
and the periastron position of the Jovian planet. 

The characteristics of each system dictated the initial 
conditions, chosen in a fine 80 by 80 grid within the hab- 
itable zone, from which the orbits were computed using 
a robust numerical method (Lie-series integration) for 1 
million years. The grid of the initial conditions of the fic- 
titious terrestrial planets was chosen to cover the whole 
habitable zone of the system, and also to model the un- 
certainties in the elements of the observed planet (s). The 
stability of orbits was assessed with two methods, namely 
the computation of the Renyi entropy as measure of the 
chaoticity of an orbit and the determination of the maxi- 
mum eccentricity of the orbit of a fictitious planet during 
its orbital evolution of 1 million years. 

We can say that our computations for such a fine grid, 
taking into account also the essential role of the MMRs, 
lead to a deeper insight concerning the dynamics of the 
five systems which we studied. We also give the percentage 
of orbits which survived in the paper (=MT) of Menou 
& Tabachnik (j2()():-{|l where they investigated all known 
extrasolar planetary systems with respect to possible ad- 
ditional terrestrial planets. We note that a direct com- 
parison of MT with the percentages of 'our' survivors 
is not useful here because of the different approaches 
used; we have emphasized the role of the MMRs and ne- 
glected possible inclinations. However, we know that ter- 
restrial planets will form within a protoplanetary disk 
thus staying with small orbital inclinations (Richardson 
et al. 120001 Lissauer ll993|l : additionally in a recent publi- 
cation (Dvorak et al. I2003afl it was shown that the incli- 
nations of the fictitious planets up to 15 degrees do not 
change the stability of orbits in the HZ. The results for 
possible stable orbits may be summarized as follows: 

— In the system Gl 777 A the stability zone for the 
motion of terrestrial planets is well inside of the HZ 
and suggests any planets residing there will survive for 
a sufficiently long time (in MT, 86.8 % of the orbits 
were found to be stable). 

— In the system 47 Uma we can also say that there is 
a good chance for planets to move inside the HZ with 



small eccentricities between the main resonances; this 
is a result which is consistent with others but not with 
MT where only 28 % remained. 

— HD 72659 turned out to be a very good candidate for 
hosting planets in the HZ; again this does not confirm 
the results of MT - they found that only 40.2 % of the 
orbits were stable. 

— The results of the computations for Gl 614 show that 
it is very unlikely that there is an additional planet 
moving in the HZ - these results are consistent with 
MT (9.2 % stable orbits) 

— For HD 4208 we can say that in the HZ there is 
enough room left for terrestrial planets and that they 
could survive for sufficiently long time; these results 
are more or less consistent with those of MT where 
50.2 % of the orbits were stable. 

New observational possibilities provided by missions 
like COROT, DARWIN or the TPF make the first 
search for terrestrial exoplanets seem possible in the next 
decades. In an ESA study the goal of the missions is sum- 
marized as follows: "To detect and study Earth-type plan- 
ets and characterize them as possible abodes of life" . In 
this sense dynamical studies like the one we present here 
should help to define promising targets for observations. 
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Appendix A: Recurrence Plots 

A.l. Introduction 

In his seminal paper Henri Poincare (|1890|l introduced the con- 
cept of recurrences in phase space, when he studied the stability 
of the solar system. Recurrence Plots (RPs) were introduced to 
visualize the recurrences of trajectories of dynamical systems 
in phase space (Eckmann et al. ll987t . These plots have proved 
to be rather useful in the analysis of time series, as they give a 
first impression about the behavior of the system under study. 
However, in order to go beyond mere visual inspection, which 
depends on the person doing the analysis and hence is always 
to some extent subjective, measures that quantify the struc- 
tures found in RPs were introduced (Webber & Zbilut lT^Mt . 
This quantification of RPs has found numerous applications in 
many fields, such as Geology (Marwan et al. 12002^ . Physiology 
(Marwan et al. 120021 Zbilut et al. l20U2t . Chmatology (Marwan 
& Kurths et al. I2UU2( . etc. Especially, it has been shown that 
it is possible to estimate dynamical invariants based on RPs, 
like the Renyi entropy and the correlation dimension (Thiel 
et al. I2003II . The estimation of these invariants by means of 
RPs has some advantages with respect to other algorithms 
usually applied to the analysis of time series. The method of 
RPs avoids, for example, the problem of the choice of embed- 
ding parameters, and it can also be applied to non-stationary 
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data (Thiel et al. I2U04(| . We concentrate on the estimation of 
the Renyi entropy, as it can be interpreted as the inverse of the 
prediction time (or Lyapunov timescale), i.e., it can be used to 
measure the prediction horizon of the system. 

A.2. Definition of Recurrence Plots and Examples 

RPs were introduced to simply visualize the behavior of tra- 
jectories in phase space (Eckmann et al. 119871 . Suppose we 
have a dynamical system represented by the trajectory {xi} 
for i = 1, . . . , A'^ in a d-dimensional phase space. Then we com- 
pute the recurrence matrix 

R^,j ^e{e-\\x,-xj\\), i,j = l...N, (A.l) 

where e is a predefined threshold and &{■) is the Heaviside 
function^. The graphical representation of Ri,j, called a 
Recurrence Plot, is obtained by encoding the value one for 
a "black" point and zero for a "white" point. The recurrence 
rate RR is defined as the percentage of black points in the RP, 
i.e. 

1 ^ 

Figure urn a) shows the RP of a sine function, i.e a circle 
in phase space. The plot is characterized by non-interrupted 
diagonal lines. Figure IA. II bl is the RP of the Rossler system 
in a chaotic regime. In this case the predominant structure 
are diagonal lines, which are interrupted. Figure lA.ll c) rep- 
resents the RP of white noise. It is homogeneous with mainly 
single points, indicating a stochastic system, as the state at 
time t + 1 is independent of the one at t. From these plots, 
we see that there is a certain connection between the length 
of diagonal lines and the ratio of determinism or predictability 
inherent to the system. This connection can be explained as 
follows: suppose that the states at t = i and t = j are neigh- 
boring, i.e. Ri^j = 1. If the system behaves predictably, then 
similar situations lead to a similar future, i.e. the probability 
for Ri+ij+i = 1 is high. For perfectly predictable systems, this 
leads to infinitely long diagonal lines (like in the RP of the sine 
function). In contrast, if the system is stochastic, the proba- 
bility for Ri+i.j+i = 1 is small, i.e., we find single points or 
short lines. If the system is chaotic, initially neighboring states 
diverge exponentially. The faster the divergence, i.e. the higher 
the Lyapunov exponent, the shorter the diagonals. In the next 
section we develop this observation and explain formally the 
relationship between the length of the diagonal lines in the RP 
and the predictability of the system. 

A. 3. Estimation of the Renyi entropy based on the RP 

In this section we recall first the definition of the Renyi entropy 
of second order and then present the mathematical relation 
between this information measure and the RPs. Let us consider 
a trajectory x{t) in a bounded d-dimensional phase space and 
the state of the system is measured at time intervals t. Let 
{1, 2, M{e)} be a partition of the attractor in boxes of size 
e. Then p{ii, ...jii) denotes the joint probability that x(t = r) 
is in the box ii, x(t — 2r) is in the box 12, and x{t = 

^ The norm used in Eq. lA.ll is in principle arbitrary. For 
theoretical reasons, it is preferable to use the maximum norm. 



It) is in the box ii. The order-2 Renyi entropy (Renvi 119701 
Grassberger 1198,^ is then defined as 

A'2 = - lim lim lim 7!- In p'^iii, . . . (A.2) 

Roughly speaking, this measure is directly related to the num- 
ber of possible trajectories that the system can take for I time 
steps in the future. If the system is perfectly deterministic, in 
the classical sense, there is only one possibility for the trajec- 
tory to evolve and hence 7^2 = 0. In contrast, one can easily 
show that for purely stochastic systems, the number of possible 
future trajectories increases to infinity so fast, that K2 — > 00. 
Chaotic systems are characterized by a finite value of K2, as 
they belong to a category between pure deterministic and pure 
stochastic systems. Also in this case the number of possible tra- 
jectories diverges but not as fast as in the stochastic case. The 
inverse of K2 has units of time and can be interpreted as the 
mean prediction horizon/time of the system. 

In the paper of Thiel et al. II2UU3II it was shown that there 
exist a direct relationship between K2 and RPs, and it is as 
follows 

lnP,"(0 ~ e-"^ exp(-A'2(e)r), (A.3) 

where P^il) is the cumulative distribution of diagonal lines in 
the RP, i.e. the probability of finding a diagonal in the RP 
of at least length I ®. D2 is the correlation dimension of the 
system under consideration (Grassberger & Procaccia 1983). 
Therefore, if we represent P^{1) in a logarithmic scale versus I 
we should obtain a straight line with slope —K2{£)t for large 
I's, which is independent of e for a rather large range in e. This 
is shown in the left panel of figure lX!2l for the chaotic Bernoulli 
map x„+i = 2xnmod{V) . One obtains K2 ~ 0.6733 for 10000 
data points, in good accordance with the values found in the 
literature (Ott ll993t . 

Finally one represents the slope of the curves for large / in 
dependence on e, where a plateau is found for chaotic systems. 
The value of this plateau determines K2 (Figure IA.2I right 
panel). If the system is not chaotic, one has to consider the 
value of the slope for a sufficiently small value of e. 

A. 4. Automatization of the Algorithm 

To compute the stability diagrams presented in this paper, 
which consist of about 13000 grid points, the estimation of K2 
had to be automated. The crucial point in the automatiza- 
tion is the estimation of the scaling region of In P^{1) vs. I and 
the plateau in 7^2 (e) vs. e. In both cases we applied a clus- 
ter dissection algorithm (Spaeth ll992t . The algorithm divides 
the set of points into distinct clusters and performs a linear 
regression in each cluster. Then the sum of all square residu- 
als is minimized in order to determine the scaling region and 
the plateau. To find both regions automatically, we used the 
following parameters: 

— We considered only diagonal lines up to length Imax = 400. 
Longer lines were excluded because of finite size effects. 

— We considered only values of Ps{l) with P^{1) > 500 for 
the same reason as in the last item. 

® Fractions of longer lines are also counted: e.g. a line of 
length 5 counts as 5 lines of length 1, 4 lines of length 2, and 
so on. 
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Fig. A.l. Prototypical examples of RPs: a) RP of a sine function, b) RP of the Rossler system in chaotic regime, c) 
RP of white noise. 




Fig. A. 2. Left panel: number of diagonal lines of at least length I versus I in the RP of the Bernoulli map for different 
values of the threshold e. The mean slope of the curves is equal to 0.6917. Right panel: Renyi entropy vs. e for the 
Bernoulli map. 



— We used 40 different values for e, corresponding to 40 
equally spaced recurrence rates RR between 1% and 95%, 
to have a well defined plateau in K2{e) vs. e. 

— We used 10000 data points of each simulated trajectory. 
The more data points one uses, the more pronounced the 
scaling regions. Note that the computation time increases 
approximately with A^^. 

— For the applied cluster dissection algorithm we had to spec- 
ify the number of clusters in each run. For the detection 
of the scaling region in InP^(Z) vs. I, we chose 2 different 
clusters and used the slope of the largest cluster. For the 
detection of the plateau in K2{e) vs. e, we chose 3 clus- 
ters and used the value of the cluster with the minimum 
absolute slope. 

These choices have proven to be the most appropriate ones for 
the estimation of the scaling regions. All these parameters are 
defaults of the computer program. The only needed input is 
the file with the trajectory data. 

Appendix B: The Lie-series IVIethod 

B.l. General Properties of the Lie-series 
Grobner 1198711 defined the Lie-operator D as follows: 

= A + (B.l) 

D is a linear differential operator; the point z — 
{zi, Z2, ■ ■ ■ , z„) lies in the n-dimensional z-space, the functions 
9i{z) are holomorphic within a certain domain G, e.g. they 



can be expanded in a converging power series. Let the function 
f{z) be holomorphic in the same region as 6i{z). Then D can 
be applied to f{z): 

Df = e,{z)^ + e^iz)^ + . . . + e^[z)^ (b.2) 

OZl_ OZ2 OZn 

If we proceed applying D to / we get 
D^f = D{Df) 

The Lie-series will be defined in the following way; 

i/=0 

Because we can write the Taylor-expansion of the exponential 
function 

e*-" = l-htZ3^ + ^D' + ^Z3' + ... (B.3) 

L{z,t) can be written in the symbolic form 
L(z,t)=e*°/(z) (B.4) 

The convergence proof of L{z, t) is given in detail in 
Grobner (ITMfl . 

One of the most useful properties of Lie-series is the 
Exchange Theorem: 
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Let F{z) be a holomorphic function in the neighborhood of 
{zi, Z2, ■ ■ ■ , Zn) where the corresponding power series expansion 
converges at the point {Zi,Z2, ■ ■ ■ , Zn); then we have: 



F{Z) = -. D'-FiZ) 
or 



F{e"')z = e*°F(2) 



(B.5) 



(B.6) 



Making use of it we can demonstrate how Lie-series solve DEs. 
Let us give the system of DEs: 



dzi 



e^{z) 



(B.7) 



with {zi, Z2, . . . , Zn). The solution of IB. 71 can be written as 



Zt = e 



(B.; 



where the are the initial conditions Zi{t — 0) and D is 
the Lie-operator as defined in IB. II In order to prove IB. 81 we 
differentiate it with respect to the time t: 



d-Zi tD f tD j-^f 

— =De ^i = e D^i 
Because of 



(B.9) 



(B.IO) 



we obtain the following result which turns out to be the 
original DE IbTtI 



dt 



^{^^) = e,{e"'^,) = e,{z,) 



(B.ll) 



D^^ ^ Dri^ = 62 

27-1 if- 



2 



D^^ = aDS, = a 7] 



For the Lie-terms even respectively odd in the order conse- 
quently one can find 

D (-1) a 5 

which leads to the solution for 2: 

z = ^ + Trj- —a ^ - —a 77 -I- —a ^ • ■ • 

Finally we get after the factorization of ^ and of 77 



which is exactly the known solution of the harmonic oscillator: 

z(t) =5 cos QfT H sin ar 

a 



B.2. A simple example 

To demonstrate the principle of the Lie-integration, we will 
show how one proceeds in the simple case of the harmonic 
oscillator, a 2"'' order DE where the solution is known: 



d^x , 2 



dt^ 



+ a X = 



(B.12) 



The first step consists in separating into two order DEs 
such that 

^^y^e.{x,y) 
^ = -a^x = 92{x,y) 

with the initial conditions z{t = 0) = ^ and y{t = 0) = 77. With 
this notation we find the Lie-operator of the form 

The solution can now be written as a Lie-series 

x = e^-^C and y = e^-^T? (B.14) 

where t — to = t and for t — this are the initial conditions. 
Being aware of the symbolic development of e^^ we can com- 
pute the first terms: 
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